RNA-Seq: differential expression analysis of controls

Libraries required

library(limma)
library(edgeR)
library(plgINS)
library(sva)
library(dplyr)
library(SummarizedExperiment)

Load salmon object

load("input/SC_controls_rnaseq_salmon.tds.RData")

Filtering gene counts

Function

counts_filter <- function(counts, nReads = 10, nSamplesPercent = 0.4) {
  counts2 <- counts[apply(counts, 1, FUN = function(x) {
    sum(x >= nReads) >= floor(length(x) * 0.4)
  }), ]
  return(counts2)
}

Counts from all data

counts.filt <- counts_filter(counts = salmon@tx.counts)

Differeitial analysis using limma

SVA

se <- SummarizedExperiment(assays = counts.filt)
se@colData <- DataFrame(salmon@phenoData)
sv.pl <- svacor(SE = se, form = ~ 0 + Group)
## Using variance-stabilizing transformation
## converting counts to integer mode
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
## Number of significant surrogate variables is:  4 
## Iteration (out of 5 ):1  2  3  4  5

limma fit

design <- model.matrix(~ 0 + SV1 + SV2 + SV3 + Group, data = sv.pl@colData)
colnames(design) <- gsub(pattern = "Group", replacement = "", x = colnames(design))

dds <- calcNormFactors(DGEList(counts = counts.filt))
v <- voom(dds, design = design)

contrast.matrix <- makeContrasts(PND15 - PND8, Adult - PND15, levels = design)
fit <- lmFit(v)

fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

DEA results

pnd8.pnd15 <- as.data.frame(topTable(fit2, coef = 1, number = Inf))
pnd8.pnd15 <- data.frame(Genes = rownames(pnd8.pnd15), pnd8.pnd15, stringsAsFactors = F)

pnd15.adult <- as.data.frame(topTable(fit2, coef = 2, number = Inf))
pnd15.adult <- data.frame(Genes = rownames(pnd15.adult), pnd15.adult, stringsAsFactors = F)

Venn Diagram

results <- decideTests(fit2)
vennDiagram(results)

vst Counts data and pData

dge <- DGEList(counts = counts.filt, group = salmon@phenoData$Group)
dge <- calcNormFactors(object = dge, method = "TMM")
cpm <- cpm(dge, log = T)
cpm <- removeBatchEffect(x = cpm, covariates = design[, 1:3], design = design[, 4:6])

data <- list(cpm = cpm, vstSV = assay(sv.pl), voomE = v$E, pData = salmon@phenoData)

save(data,
  file = "./output/data_pData_tx.RData", compress = T,
  compression_level = 3
)

Results

dea.list <- list(
  `PND8 vs PND15` = as.DEA(pnd8.pnd15),
  `PND15 vs Adult` = as.DEA(pnd15.adult)
)

normCounts <- v$E

voomEList <- v

dea.limma <- list(
  `PND8 vs PND15` = pnd8.pnd15,
  `PND15 vs Adult` = pnd15.adult
)

Save RData files

save(dea.list,
  file = "./output/tx_dea_SC_Controls.DEA.RData", compress = T,
  compression_level = 3
)

save(normCounts,
  file = "./output/tx_normCounts_voom_SC_Controls.RData", compress = T,
  compression_level = 3
)

save(voomEList,
  file = "./output/tx_voom_EList_SC_Controls.RData", compress = T,
  compression_level = 3
)

save(dea.limma,
  file = "./output/tx_limma_SC_Controls.RData", compress = T,
  compression_level = 3
)

PCA plots

Raw counts filtered

plPCA(x = counts.filt, samples_data = salmon@phenoData, colorBy = "Group", add.labels = FALSE)

SVA (VST)

plPCA(x = sv.pl@assays@data$corrected, samples_data = salmon@phenoData, colorBy = "Group", add.labels = FALSE)

CPM

plPCA(x = cpm, samples_data = salmon@phenoData, colorBy = "Group", add.labels = FALSE)

Voom

plPCA(x = v$E, samples_data = salmon@phenoData, colorBy = "Group", add.labels = FALSE)

SessionInfo

devtools::session_info() %>%
  details::details()

─ Session info ───────────────────────────────────────────────────────────────
 setting  value                       
 version  R version 3.6.2 (2019-12-12)
 os       Ubuntu 16.04.6 LTS          
 system   x86_64, linux-gnu           
 ui       X11                         
 language (EN)                        
 collate  en_US.UTF-8                 
 ctype    en_US.UTF-8                 
 tz       Europe/Zurich               
 date     2020-07-08                  

─ Packages ───────────────────────────────────────────────────────────────────
 package              * version    date       lib
 acepack                1.4.1      2016-10-29 [1]
 annotate               1.64.0     2019-10-29 [1]
 AnnotationDbi          1.48.0     2019-10-29 [1]
 assertthat             0.2.1      2019-03-21 [1]
 backports              1.1.7      2020-05-13 [1]
 base64enc              0.1-3      2015-07-28 [1]
 Biobase              * 2.46.0     2019-10-29 [1]
 BiocGenerics         * 0.32.0     2019-10-29 [1]
 BiocParallel         * 1.20.1     2019-12-21 [1]
 Biostrings             2.54.0     2019-10-29 [1]
 bit                    1.1-15.2   2020-02-10 [1]
 bit64                  0.9-7      2017-05-08 [1]
 bitops                 1.0-6      2013-08-17 [1]
 blob                   1.2.1      2020-01-20 [1]
 bookdown               0.18       2020-03-05 [1]
 callr                  3.4.3      2020-03-28 [1]
 checkmate              2.0.0      2020-02-06 [1]
 cli                    2.0.2      2020-02-28 [1]
 cluster                2.1.0      2019-06-19 [1]
 colorspace             1.4-1      2019-03-18 [1]
 crayon                 1.3.4      2017-09-16 [1]
 crosstalk              1.1.0.1    2020-03-13 [1]
 data.table             1.12.8     2019-12-09 [1]
 DBI                    1.1.0      2019-12-15 [1]
 DelayedArray         * 0.12.3     2020-04-09 [1]
 desc                   1.2.0      2018-05-01 [1]
 DESeq2               * 1.26.0     2019-10-29 [1]
 devtools               2.3.0      2020-04-10 [1]
 digest                 0.6.25     2020-02-23 [1]
 dplyr                * 1.0.0      2020-05-29 [1]
 edgeR                * 3.28.1     2020-02-26 [1]
 ellipsis               0.3.1      2020-05-15 [1]
 evaluate               0.14       2019-05-28 [1]
 fansi                  0.4.1      2020-01-08 [1]
 farver                 2.0.3      2020-01-16 [1]
 foreign                0.8-76     2020-03-03 [1]
 Formula                1.2-3      2018-05-03 [1]
 fs                     1.4.1      2020-04-04 [1]
 genefilter           * 1.68.0     2019-10-29 [1]
 geneplotter            1.64.0     2019-10-29 [1]
 generics               0.0.2      2018-11-29 [1]
 GenomeInfoDb         * 1.22.1     2020-03-27 [1]
 GenomeInfoDbData       1.2.2      2019-11-18 [1]
 GenomicRanges        * 1.38.0     2019-10-29 [1]
 GEOquery               2.54.1     2019-11-18 [1]
 ggplot2              * 3.3.1      2020-05-28 [1]
 glue                   1.4.1      2020-05-13 [1]
 gridExtra              2.3        2017-09-09 [1]
 gtable                 0.3.0      2019-03-25 [1]
 Hmisc                  4.4-0      2020-03-23 [1]
 hms                    0.5.3      2020-01-08 [1]
 htmlTable              1.13.3     2019-12-04 [1]
 htmltools              0.4.0      2019-10-04 [1]
 htmlwidgets            1.5.1      2019-10-08 [1]
 httr                   1.4.1      2019-08-05 [1]
 IRanges              * 2.20.2     2020-01-13 [1]
 jpeg                   0.1-8.1    2019-10-24 [1]
 jsonlite               1.6.1      2020-02-02 [1]
 knitr                  1.28       2020-02-06 [1]
 lattice                0.20-41    2020-04-02 [1]
 latticeExtra           0.6-29     2019-12-19 [1]
 lazyeval               0.2.2      2019-03-15 [1]
 lifecycle              0.2.0      2020-03-06 [1]
 limma                * 3.42.2     2020-02-03 [1]
 locfit                 1.5-9.4    2020-03-25 [1]
 magrittr               1.5        2014-11-22 [1]
 Matrix                 1.2-18     2019-11-27 [1]
 matrixStats          * 0.56.0     2020-03-13 [1]
 memoise                1.1.0.9000 2020-05-06 [1]
 mgcv                 * 1.8-31     2019-11-09 [1]
 munsell                0.5.0      2018-06-12 [1]
 nlme                 * 3.1-147    2020-04-13 [1]
 nnet                   7.3-13     2020-02-25 [1]
 pillar                 1.4.4      2020-05-05 [1]
 pkgbuild               1.0.8      2020-05-07 [1]
 pkgconfig              2.0.3      2019-09-22 [1]
 pkgload                1.1.0      2020-05-29 [1]
 plgINS               * 0.1.5      2020-07-08 [1]
 plotly               * 4.9.2.1    2020-04-04 [1]
 plyr                   1.8.6      2020-03-03 [1]
 png                    0.1-7      2013-12-03 [1]
 prettyunits            1.1.1      2020-01-24 [1]
 processx               3.4.2      2020-02-09 [1]
 ps                     1.3.3      2020-05-08 [1]
 purrr                  0.3.4      2020-04-17 [1]
 R6                     2.4.1      2019-11-12 [1]
 RColorBrewer           1.1-2      2014-12-07 [1]
 Rcpp                   1.0.4.6    2020-04-09 [1]
 RCurl                  1.98-1.2   2020-04-18 [1]
 readr                  1.3.1      2018-12-21 [1]
 remotes                2.1.1      2020-02-15 [1]
 rlang                  0.4.6      2020-05-02 [1]
 rmarkdown              2.1        2020-01-20 [1]
 rmdformats             0.4.0      2020-06-07 [1]
 rpart                  4.1-15     2019-04-12 [1]
 rprojroot              1.3-2      2018-01-03 [1]
 RSQLite                2.1.4      2019-12-04 [1]
 rstudioapi             0.11       2020-02-07 [1]
 S4Vectors            * 0.24.4     2020-04-09 [1]
 scales                 1.1.1      2020-05-11 [1]
 sessioninfo            1.1.1      2018-11-05 [1]
 SRAdb                  1.48.2     2019-12-24 [1]
 stringi                1.4.6      2020-02-17 [1]
 stringr                1.4.0      2019-02-10 [1]
 SummarizedExperiment * 1.16.1     2019-12-19 [1]
 survival               3.1-12     2020-04-10 [1]
 sva                  * 3.34.0     2019-10-29 [1]
 testthat               2.3.2      2020-03-02 [1]
 tibble                 3.0.1      2020-04-20 [1]
 tidyr                  1.1.0      2020-05-20 [1]
 tidyselect             1.1.0      2020-05-11 [1]
 usethis                1.6.1      2020-04-29 [1]
 vctrs                  0.3.1      2020-06-05 [1]
 viridisLite            0.3.0      2018-02-01 [1]
 withr                  2.2.0      2020-04-20 [1]
 xfun                   0.13       2020-04-13 [1]
 XML                    3.99-0.3   2020-01-20 [1]
 xml2                   1.3.2      2020-04-23 [1]
 xtable                 1.8-4      2019-04-21 [1]
 XVector                0.26.0     2019-10-29 [1]
 yaml                   2.2.1      2020-02-01 [1]
 zlibbioc               1.32.0     2019-10-29 [1]
 source                          
 CRAN (R 3.6.1)                  
 Bioconductor                    
 Bioconductor                    
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 Bioconductor                    
 Bioconductor                    
 Bioconductor                    
 Bioconductor                    
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 Bioconductor                    
 CRAN (R 3.6.1)                  
 Bioconductor                    
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 Bioconductor                    
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 Bioconductor                    
 Bioconductor                    
 CRAN (R 3.6.1)                  
 Bioconductor                    
 Bioconductor                    
 Bioconductor                    
 Bioconductor                    
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 Bioconductor                    
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 Bioconductor                    
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 Github (r-lib/memoise@4aefd9f)  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 local                           
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 Github (juba/rmdformats@94cd7a3)
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 Bioconductor                    
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 Bioconductor                    
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 Bioconductor                    
 CRAN (R 3.6.2)                  
 Bioconductor                    
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.2)                  
 CRAN (R 3.6.1)                  
 Bioconductor                    
 CRAN (R 3.6.2)                  
 Bioconductor                    

[1] /home/ubuntu/R/x86_64-pc-linux-gnu-library/3.6
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library